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ABSTRACT 

Global solutions of optically thick advective accretion disks around black holes are 
constructed. The solutions are obtained by solving numerically a set of ordinary differential 
qq ' equations corresponding to a steady axisymmetric geometrically thin disk. We pay special 

. attention to consistently satisfy the regularity conditions at singular points of the equations. 

For this reason we analytically expand a solution at the singular point, and use coefficients of 
the expansion in our iterative numerical procedure. We obtain consistent transonic solutions in 
a wide range of values of the viscosity parameter a and mass accretion rate. We compare two 
different form of viscosity: one takes the shear stress to be proportional to the pressure, while 
the other uses the angular velocity gradient-dependent stress. 

We find that there are two singular points in solutions corresponding to the pressure- 
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proportional shear stress. The inner singular point locates close to the last stable orbit around 
■ black hole. This point changes its type from a saddle to node depending on values of a and 

accretion rate. The outer singular point locates at larger radius and is always of a saddle-type. 
We argue that, contrary to the previous investigations, a nodal-type inner singular point does 
not introduce multiple solutions. Only one integral curve, which corresponds to the unique 
global solution, passes simultaneously the inner and outer singular points independently of the 
type of inner singular point. Solutions with the angular velocity gradient-dependent shear stress 
have one singular point which is always of a saddle-type and corresponds to the unique global 
solution. The structure of accretion disks corresponding to both viscosities are similar. 
Subject headings: Accretion, accretion disks — black hole physics — hydrodynamics 
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1. Introduction 

Accretion discs are formed when the matter with a large angular momentum is falling into a black hole 
or another gravitating body. The well known objects where the accretion disks are found are protostar 
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ncbulac, binary X-ray sources, cataclysmic variables, active galactic nuclei and others. In this paper we 
discuss accretion disks around black holes. The standard model of geometrically thin accretion disk has 
been developed by Shakura (1972), Novikov & Thorne (1973) and Shakura & Sunyaev (1973), and has 
played a significant role in the development of accretion theory (see Pringle 1981; Frank, King, & Rain 1992; 
Kato, Mineshige, & Fukue 1998 for reviews). The standard model bases on the vertically averaged approach 
to equilibrium, and a suggestion of the local thermal balance in which the viscous heating of the gas is 
balanced by the local radiative cooling. Non-local effects, like the radial advection of thermal energy and 
the transonic nature of accretion flow, are neglected in the standard model. This simplified approach allows 
to reduce the general problem to a set of algebraical equations. Such a simple description becomes possible 
due to an approximate parameterization of the viscosity stress tensor with one non-zero component, 

t rip = -aP, (1) 

suggested by Shakura (1972). The standard model gives a satisfactory appropriate solution of the problem 
at low accretion rates M <; WLsdd/c 2 , where Lsdd is the Eddington luminosity. 

Simplified solution with inclusion of the advective terms into equations described the vertically 
integrated models of accretion disks was obtained by Paczyhski & Bisnovatyi-Kogan (1981). This approach 
with some modifications have been used by many researchers to study transonic accretion flows around 
black holes (Muchotrzeb & Paczyhski 1982; Muchotrzeb 1983; Matsumoto ct al. 1984; Fukue 1987; 
Abramowicz et al. 1988; Chen & Taam 1993; Beloborodov 1998). The importance of the transonic nature 
of the accretion flows on the disk structure has been emphasized by Hoshi & Shibazaki (1977), Liang & 
Thompson (1980) and Abramowicz & Zurek (1981), and later studied in more details by Abramowicz & 
Kato (1989). 

Despite a significant progress in the study of optically thick accretion disks obtained during almost 
three decades there are a number of unsolved problems still posed in the theory. The problems are connected 
with a possible non-uniqueness of a solution at a ^ 0.01 and a non-standard behavior of a singular point 
type. It was reported by Matsumoto et al. (1984), Muchotrzeb-Czerny (1986) and Abramowicz et al. 
(1988) that in the case of viscosity prescription (1) the singular point changes its type from a saddle to 
node when one increases a. The presence of the nodal- type singular point leads to creating of a possibility 
of multiple solutions as the authors have claimed. A similar change of the singular point type was reported 
by Chen & Taam (1993), who used the angular velocity gradient-dependent viscous stress, 

t r4> = pvr—, (2) 

where v is the kinematic viscosity coefficient defined by (14). Narayan, Kato, & Honma (1997) have 
compared two forms of viscosity (1) and (2) in the case of radiatively inefficient advection-dominated 
accretion flows. They concluded that the structures of flows corresponded to both viscosities are similar at 
a < 0.15. 

In this paper we show that the mentioned problems have been created by several inconsistencies in 
the preceding studies. Some problems are connected with an inaccurate averaging of the equations over 
a disk thickness (Chen & Taam 1993), another ones appear due to an incomplete investigation of the 
singular points (Abramowicz et al. 1988). We have found that in the case of viscosity prescription (1) 
a set of equations describing the vertically averaged advective accretion disks has two singular points, 
independed of a and accretion rate. Note, that multiplicity of singular points in solutions for accretion 
flows in Paczyhski- Wiita potential (3) was revealed by Fukue (1987), Chakrabarti & Molteni (1993) and 
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Chakrabarti (1996) in a somewhat different context. We have shown, that at a <J 0.01 the inner and 
outer (with respect to the black hole location) singular points are of the saddle type, and only one integral 
curve ( "separatrix" ) which crosses the inner point simultaneously crosses the outer one. This separatrix 
corresponds to the unique global solution which is determined by two parameters, a and m = Mc 2 /L Edd , 
for a given black hole mass. In Figure la the structure of integral curves is schematically represented in the 
vicinity of the global solution which is shown by the thick line. At larger a > 0.1 the inner singular point 
is changed to a nodal-type one, while the outer point remains of a saddle-type. There is still one integral 
curve which goes continuously through both singular points providing a unique global solution, as it is 
shown in Figure lb. 

In the case of viscosity prescription (2) we have found that there is only one singular point which is 
always a saddle, and only one physical solution which passes through this point exists. Solutions which 
correspond to both forms of viscosity (1) and (2) are very close at low a limit, a <; 0.1. 

We have developed a numerical method to solve the set of equations describing the vertically averaged 
advective accretion disks. The method is based on the standard relaxation technique and explicitly uses 
conditions at the inner singular point and its vicinity. We have obtained these conditions by expanding a 
solution into power series around the singular point. Such a modification of the method allows to construct 
solutions which smoothly pass the singular points and satisfy the regularity conditions at these points with 
high computer precision in wide range of parameters a and m. 

The paper is organized as follows. In §2 we formulate a mathematical approach to the problem, write 
a set of equations, and formulate boundary conditions. In §3 we investigate critical points and discuss 
behavior of physical values in their vicinity. In §4 we describe our numerical results and discuss them in §5. 
Details of the numerical method and explicit expansion of physical quantities in the vicinity of the critical 
points are represented in Appendixes A and B, respectively. 



2. Problem formulation 

We consider a steady geometrically thin accretion disk around a non-rotating black hole. For simplicity, 
we use the pseudo-Newtonian approach to describe the disk structure in the vicinity of a black hole. In 
the approach the general relativistic effects are simulated by using Paczyhski-Wiita potential (Paczyhski & 
Wiita 1980) 

where M is the black hole mass and 2r g = 2GM/c 2 is the gravitational radius. The disk self-gravity is 
neglected. 

A general problem of investigation of two-dimensional structure of the accretion disks (in the radial 
and vertical directions) can be reduced to a one-dimension problem by averaging the disk structure in the 
vertical direction. In this formulation equations which are described the radial disk structure are written 
for the midplane density p, pressure P, radial velocity v and angular velocity fi. The mass conservation 
equation takes the form, 

M = inrhpv, (4) 
where M is the accretion rate, M > 0, and h is the disk half-thickness, which is expressed in terms of the 
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isothermal sound speed c s = \fPj~p of gas, 

The equations of motion in the radial and azimuthal directions are 

dv _ ldP 
dr p dr 



1 dP \ 
w— = — + (0 - f2|r)r, 



T-T" + -r( r ht r<p) = 0, 
47rdr dr v 



(5) 

(6) 
(7) 



where £Ik is the Keplerian angular velocity, Vl 2 K = GM/r(r — 2r s ) 2 , £ = fir 2 is the specific angular 
momentum and t r $ is the (r, </>)-component of the viscous stress tensor. Other components of the stress 
tensor are assumed to be negligiblly small. 



The vertically averaged energy conservation equation can be written in the formjj], 



where 



F„ 



r adv . . 

Zirr 
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F~ 



dE d ( 1 
dr dr \p 



F = ht r(j> r—, 
dr 

2aT 4 c 



3nph 



(8) 

(9) 
(10) 

(11) 



are the advective energy flux, the viscous dissipation rate and the cooling rate per unit surface, respectively, 
T is the midplane temperature, k is the opacity and a is the radiation constant. 



The equation of state for accretion matter consisted of a gas-radiation mixture is 

3 p 

where 1Z is the gas constant. The specific energy of the mixture is 



E = -TIT 
2 



aT 4 
P 



(12) 



(13) 



We consider two prescriptions of viscosity in our models. In one case we adopt a simple relation (1) 
between the viscous stress and pressure. In another case we assume the angular velocity gradient-dependent 
viscous stress (2), where the viscosity v is taken in the form 

2 

v = -ac s h. (14) 
In the limit O — ► VLk both prescriptions (1) and (2) coincide. 



The vertical averaging in equation (8) have been done in different way by different authors (compare e.g. Shakura & 
Sunyaev 1973, and Abramowicz et al. 1988). Our choice of coefficients in (9)-(ll), following Chen & Taam (1993) may be 
not the optimal one. Aposteriory analysis had shown that using 4 in the denominator of (9) instead of 2, would be a more 
consistent choice, but this change has a little influence on our numerical results. 
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Integrating equation (7) we obtain 

r 2 ht r<p = - — (e-£ in ), (15) 

where £{ n is an integration constant and has a meaning of the specific angular momentum of accreting 
matter at the black hole horizon. Depending on the used viscosity prescription (1) or (2) expression (15) 
results in an algebraical equation or a first order differential equation, respectively. So, in the case of 
viscosity prescription (1) we have only two first order differential equations (6) and (8), which require to 
set two parameters as boundary conditions. In the case of viscosity prescription (2) we obtain additional 
differential equation from (15) and have to set three parameters as boundary conditions. The integration 
constant £i„ is chosen to obtain a global transonic solution with a subsonic part at large radii and a 
supersonic part close to the black hole horizon. 



3. Investigation of singular points 

3.1. aP viscosity prescription 

We consider first the case of viscosity prescription (1). From (15) we obtain the algebraical expression 
for fl, 

Q=%+c^. (16) 
r vr 

Using (16) the system of equations (6) and (8) can be reduced to the following form, 

r^ = {l- M 2 ) — L + l- r -JL + 5— ^r 2 , (18) 

c s Di il K cj 



where 



Q K c 2 J \ 2'' 4-3/3 M 2 J V n K J \ 2^4-3/3 

vr 2 A4 Z m c s r g 

In equations (17)-(20) we use the following notations: v' = dv/dr, c' s = dc s /dr, Q' K = dflx/dr, [3 = 1ZT/c 2 
and M — v/c s . From (4) and (16) it follows the algebraical equation for (3, 

^_ (1 _ ffl f*Gs*_o. <2i) 

47r arvc' s 

The equation D\ = defines singular points of the differential equations (17) and (18), and can be 
reduced to the following form, 
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Equation (22) is a quadratic equation with respect to M 2 and has two positive roots, which correspond to 
two singular points: 




a 2 + 8 _ ± W ( a 2 + 8 - -^lr) -2a? (7- 1 " ' 



4-3/3, VV 4-3/3J V 2"°4-3# 



r-?A 1+A 



2^4 -3ft 



-1 



(23) 



where /3 S is the value of [3 taken at the singular point. 
In a <C 1 limit we have: 



^-1.-^(7-1^1 



and 



• m =-t( 8 -j^i 



The first singular point, in which .M s = .Mi, locates close to the black hole last stable orbit at r — 6r g . The 
corresponding values of Mi are 1.118 and 1.069 for the gas pressure supported (/? = 1) and the radiation 
pressure supported {(3 = 0) accretion flows, respectively. This point is an analogy of the singular point 
in a spherical flow, where the point divides the subsonic and supersonic regions of accretion flow. The 
second singular point, in which M s = M2, located at larger radius, is the result of simplified viscosity 
prescription (1). We will use the notations (r s )j„ and (r s ) out for positions of the inner and outer singular 
points, respectively. 

At the singular points the numerator N\ and denominator Di must simultaneously vanish to provide a 
regular behavior for a global solution. The type of the singular points must be consistent with a transonic 
nature of solution. For example, a spiral-type singular point does not satisfy the latter requirement, but a 
saddle or nodal- type point does it. Detailed analysis of topology in vicinity of singular points was done for 
thin accretion disks under isothermal approximation by Abramowicz & Kato (1989). They showed that 
saddle, nodal or spiral types are formally possible, but only saddle and nodal points are physically relevant. 
This study had confirmed the previously obtained numerical results by Matsumoto et al. (1984). The latter 
authors demonstrated in a framework of the isothermal accretion disks that the type of the inner singular 
point is defined by value of l in . At larger l in the point is a spiral, at smaller l in there is no inner singular 
point at all, and only unique choice of £i n of moderate values corresponds to a saddle or nodal-type singular 
point. The choice between saddle or nodal- type singular points can be done only by constructing a global 
model of the disk. 



3.2. Sl-gradient-dependent viscous stress 

In the case of viscosity prescription (2) the differential equations (6), (8) and (15) can be reduced to 
the following form, 

n' 3 n K rv ( t m \ 

r 7T = -o— -sT U-oIa > ( 26 ) 



Q 2 ac 2 s V Mr 2 , 
v' N 2 

r< = (1 - M',^ + 1 -r£k + £l_* r2> (28) 

C s D 2 ilK cj 
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wherc fi' = dQ/dx, (3 is defined by (21) and 



n' K n 2 -n 2 K2 \f^ 3 a 1 + /? \ , 3 n 2 n K r 3 v £ in \ 2 



tt' K \ f, , 3 1 -/? \ 1 -/3ft K r 2 



ft _ (A ,. 1) ( 7 _l^).( 1+ l,w ) . 

There is only one singular point of equations (26)-(28) defined by the equation D2 = 0. The point is 
an analogy to the inner singular point discussed in §3.1. To be consistent with §3.1 we use notation (r s )i n 
for the position of the point. At (r s )i n we have 

Note, that the expression for A4\ given by (24) coincides with (31). The latter could mean that the 
properties of global solutions in the case of viscosity prescriptions (1) and (2) are very similar in the inner 
part of flow at the limit of small viscosity a<l. Our numerical results confirm this conclusion. 

Abramowicz & Kato (1989) studied analytically the type of singular point in the isothermal disks in 
the case of viscosity prescription (2) . They showed that the point is always a saddle, and there is no case of 
a node. This conclusion differs from one obtained in the case of viscosity prescription (1). Our numerical 
models confirm this dependence of the singular point type on a form of viscosity. 



4. Numerical results 

To be used in the numerical method the sets of differential equations (17)-(18) and (26)-(28) have 
been re- written in the dimensionless form using the following dimensionless quantities: f = r/r g , v = v/c, 
c s = c s /c, 17 = Qr g /c, £ = £/r g c, m = M/M Q . In the subsequent discussions we will use these dimensionless 
quantities skipping in the notations the 'tilde' mark. The used method is described in Appendix A. We 
have calculated a number of models varied by the viscosity prescriptions and parameters a and tin. The 
black hole mass m contributes into the dimensionless equations in the combinations 

clc^kGMq m 

The parameter D was taken to be inversely proportional to tin with m = 10, 1Z = 1.65 • 10 8 erg g~ x K~ x and 
k = 0.4 cm 2 g~ 1 . The numerical grid covers the radial range form r in located at the inner singular point, 

Tin = ( r s)in, till T " ou t ~ 10 Tg. 

We discuss first the influence of the numerical outer boundary conditions on our models. We have 
found that the models are insensitive to the specific values of the outer boundary conditions. By fixing 
a and rh the unique global transonic solution is fully determined. This solution also uniquely determines 
the outer boundary values: two or three values depending on the used viscosity prescription (1) or (2), 
respectively. In general, we do not know a priori these 'correct' boundary values which the global solution 
passes through, and consequently, our assumed numerical boundary values are quite arbitrary. But, this 
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values should be close enough to the 'correct' ones due to reason of numerical stability. Also, the 'correct' 
outer values are very close, but not exactly equal, to the values obtained from the standard model for the 
Keplerian accretion disks. Calculations show that all numerical solutions which have the same a and to, but 
different boundary values at different r out , converge to the 'common' solution which is not affected by the 
outer boundary. This 'common' solution represents the global solution which we seek. Significant differences 
between some numerical solution and the 'common' one (with relative errors ;> 10~ 4 ) are observed only in 
2 — 3 grid points before the last outer point at r out . Such a behaviour of the numerical solutions can be 
explained by special properties of difference equations (Al) at £ ~ 1. At small e the method is unstable. 

Each model is characterized by value of [see eq.(15)] which has a sense of a specific angular moment 
of matter infalling into black hole. Figure 2 [panels (a) and (b)] shows the dependence of £i n on accretion 
rate to for three values of a = 0.01, 0.1 and 0.5, and two forms of viscosity prescription (1) and (2). At low 
to <; 1 the value of £ in is independent of to, but weakly varies with a. In the low a case, a = 0.01 and 0.1, 
the values of £ in are close to the minimum value of the Keplerian angular momentum, (£ K ) min = 3.6742. At 
high to ;> 0.1 the values of £i n deviate from [Ik)™,™ to larger or smaller values depending on a. In the case 
a = 0.01 and 0.1 one can see only minor differences between models with different forms of viscosity. But, 
for large a = 0.5 the difference in values of £i n increases. Unfortunately, we had been able to calculated only 
a limited number of models in the case of viscosity prescription (2) due to technical reason (see Appendix A 
for details), and our comparison of both prescriptions is not complete in this respect. 

Figure 2 [panels (c) and (d)] shows locations of the inner singular points (r s )i n as a function of to for 
different values of a and two different viscosity prescriptions. Similar to the case of £i n the models at low 
to show a weak dependence of {r s ) in on m. In the low a models (squares and circles in Figure 2) the values 
of {r s )in are close to the location of the black hole last stable orbit at r = 6. At high to ;> 0.1 the values of 
{r s )in are decreasing functions of to in the case of low a — 0.01 and 0.1, and non-monotonically behave in 
the case of a = 0.5 (triangles in Figure 2). 

Figure 3 shows locations of the outer singular points {r s ) out in the case of viscosity prescription (1) 
as a function of to for different values of a. Values of (r s ) out are an increasing function of m and show a 
power-law behaviour at to ^ 3. It is interesting to note that values of (r s ) out are almost independent of a. 

Examples of the specific angular momentum distribution £{r) are shown in Figure 4 for to = 160 
and three values of a — 0.01 (short-dashed line), 0.1 (solid line) and 0.5 (dotted line). The distributions 
correspond to viscosity prescription (1). The location of the inner singular points are indicated by the 
correspondent points on the curves. The Keplerian angular momentum is displayed by the long-dashed line 
for comparison. Only the low viscosity model with a = 0.01 has a super-Kcplerian part in £(r). Models 
with larger viscosity are everywhere sub-Keplerian. Note, that the singular point in the low viscosity model 
(short-dashed line) locates in the inner sub-Keplerian part of the disk. 

Figure 5 shows the dependence of (3 S on to at the inner singular points. The change of value of (3 form 
1 to corresponds to the change of a state from the gas pressure to radiative pressure dominated one. The 
thin disks with (3 ~ 1 are locally stable, whereas the parts of the disk in which (3 ~ are thermally and 
viscously unstable at to ^ 100 (Pringle, Rees, & Pacholczyk 1973). At larger to ;> 100 the instability can be 
suppressed by the advection effect (Abramowicz et al. 1988). We have found a weak dependence of (3 s (m) 
on the assumed viscosity prescriptions. 

Using analysis discussed in Appendix B we have determined a type of singular points in our numerical 
solutions. In Figures 2-5 the saddle-type points are indicated by the solid squares, circles and triangles. 
The nodal-type points are represented by the corresponding empty dots in the same figures. In the case of 
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viscosity prescription (1) the solutions have two singular points (see §3.1). The inner singular points, (r s ) in , 
can be saddles or nodes depending on values of a and m. Note that the change of type from a saddle to 
nodal one does not introduce any features in solutions. The outer singular points, {r s ) out , are always of a 
saddle-type. In the case of viscosity prescription (2) the solutions have only inner singular points (see §3.2) 
which are always of a saddle-type. 

In models with low m <J 16 and low a <J 0.1 the values of r s and ii n are quantitatively very close to 
the last stable orbit location (rj„ = 6) and value of £i n = {(K)min] assumed in the standard model 
(Shakura & Sunyaev 1973). The radial structure in our models is also very close to the one for the standard 
model in the same range of m and a. Such a good coincidence means that the advective terms in equations 
(6) and (8) are negligiblly small in considered models. However, the high a models show quite significant 
deviation from the standard model independently of m (see Figure 2). 

At high accretion rates, m ;> 16, the effect of advection becomes significant. We illustrate it by 
calculating the luminosity L of disk in the case of viscosity prescription (1), 



where F~ is given by (11). Figure 6 shows calculated dependences of L/LEdd on m for different values 
of a — 0.01, 0.1, 0.5 (short-dashed, solid and dotted lines, respectively). There is a simple linear relation 
L/LEdd = V 1 ™ m t ne standard model, in which the advection is neglected. The radiative efficiency 77 is a 
constant and equals 77 = 1/16 in the case of gravitational potential (3). We plot this relation by the straight 
long-dashed line in Figure 6. One can clearly see from the figure that effect of advection results in reduction 
of luminosities with respect to the one for the standard model at m ^ 16. There is a weak dependence of 
the luminosity on value of viscosity in disks. 

Finally note, that our numerical solutions corresponding to viscosity prescription (1) have some 
resemblance to the results of Abramowicz et al. (1988), but they show important quantitative differences, 
especially for large a and m. 



We have obtained unique solutions for structure of advective accretion disk in a wide range of accretion 
rates and a-parameters. Both viscosity prescriptions (1) and (2) have been investigated. The solutions 
corresponding to both prescriptions are very close for a <J 0.1, and begin to differ at larger a. This is 
connected, probably with larger deviation of the angular velocity fl from the Keplerian one, SIk , leading to 
larger difference between t r< f, in both prescriptions. Unfortunately, our comparison of the prescriptions is 
not complete due to technical problems in calculation of the high viscosity models in the case of viscosity 
prescription (2). 

The main difference of the present study from the previous ones is in using more sophisticated numerical 
technique which accurately treats the regularity conditions in the inner singular point of equations. We have 
performed an analytical expansion at the singular point to calculate the derivatives of physical quantities. 
These derivatives help us to find the proper integral curve passing through the singular point. The approach 
allow us to avoid numerical instabilities and inaccuracies, appearing when only variables at the singular 
point, but not its derivatives, are included into a numerical scheme. 




(32) 



5. 



Discussion 



We have found different behaviour of integral curves depending on used viscosity prescription. In the 
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case of viscosity prescription (1) there are two singular points located at (r s )i n and (r s ) out . The inner point, 
(r s )in, locates close to the last stable black hole orbit (see Figure 2), and is an analogy of the singular 
point in spherical flow, where the point divides the subsonic and supersonic regions. The location of the 
outer point, (r s ) out , is determined by the accretion rate (see Figure 3). At low a < 0.1 both points are 
of a saddle-type. Only one integral curve ( "separatrix" ) simultaneously crosses two saddle- type points, as 
it is shown in Figure la, and corresponds to the global solution which smoothly connects the supersonic 
innermost region of the accretion disk and the subsonic outer (formally at r = oo) parts. For larger a ^ 0.1 
the inner singular point changes its type to a node. There was suggestion by Muchotrzeb-Czerny (1986) 
and Abramowicz et al. (1988) that there is no unique solution in this case, because of all integral curves 
cross the node. Existence of a unique separatrix crossed simultaneously both singular points preserves a 
uniqueness of the solution in this case (see Figure lb). The conclusion of Muchotrzeb-Czerny (1986) and 
Abramowicz et al. (1988) is probably connected with their neglection of outer singular points inherited to 
the problem. 

Matsumoto et al. (1984) used slightly different form of equations (6) and (8), and they found that in 
this case only one singular point exists and changes type from a saddle to node. The difference with respect 
to our results arises because of using different form of the pressure gradient force [the first term on the right 
hand side of equation (6)]. We use the pressure and density taken at the equatorial plane in this term, 
whereas Matsumoto et al. (1984) used the vertically averaged quantities in it. In the latter case the term 
has the following form, 



The vertically integrated approach (33) introduces the difference because in this case the free terms with 
a 2 in (20) are absent. Formally, it corresponds to location of the outer singular point at infinity, where 
M2 = 0. Such a visible difference is not qualitatively important for the physical solution, because conditions 
at the outer singular point are only shifted to infinity, and the integral curve itself has little changes. Thus, 
similar to our results, in the approach by Matsumoto et al. (1984) the inner critical points of both types, a 
saddle or node, correspond to a unique solution. 

In the case of viscosity prescription (2) we have found that one singular point exist. The point is 
always of a saddle-type and determines a unique solution. Another results had been obtaining by Chen & 
Taam (1993). They also found that equations have one singular point, but the point changes its type from a 
saddle to nodal one, depending on a and the accretion rate. It is not clear why such a result was obtained. 
There are two main differences in our equations (6) and (8), and those used by Chen & Taam (1993). First, 
they used the same vertical averaging for the equation of motion as Matsumoto et al. (1984). Second, they 
used the vertically averaged energy equation which corresponds to the inappropriate polytropic relation, 
? oc S 7 , when one neglects terms corresponding to the viscous heating and radiative cooling. Here 7 
is the effective adiabatic index, and other notations are similar to those used in (33). Our equation (8) 
corresponds to the correct polytropic relation, F a p 7 . It could be that one of the mentioned differences 
results in the change of critical point type. 
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Y,dr 



where 




and 




(33) 



A. Numerical method 



We use the finite-difference method to solve the systems of ordinary differential equations discussed in 
§3. The method has some resemblance to that used by Igumenshchev, Abramowicz & Novikov (1998). In 
this approach the problem is reduced to a solution of a system of non-linear algebraical equations written 
for a each pair of neighboring numerical grid points. The numerical grid {r^} extends over about three 
orders of magnitude in the radial direction. We look for a numerical solution in which the location of the 
inner grid point n coincides with the location of the singular point (r s )in near black hole. In our method 
we approximate the differential equation dy/dr = f(r,y) by the following finite differences, 

Vi-Vi-i =e/ ._ 1 _ ( i_ g)/ z = l,2,...,I, (Al) 

where the function y(r) must be replaced by v(r), c^(r), and additionally by f2(r) in the case of viscous 
prescription (2), e is a parameter, I is the number of grid points, and the lower indices indicate the 
corresponding grid point. The value of parameter e is chosen to provide a stability of the numerical scheme. 
We use e = 1 in most of the cases. 

We use the Newton- Raphson iteration scheme to solve the set of equations (Al). In general formulation 
equations (Al) can be represented by K functional relations, involving K variables yu, 

F k (y 1 ,y 2 ,...,y K ) = 0, k = l,2,...,K, (A2) 

or in the vector notation, F(y) = 0. The (n + l)-iteration improvement of an approximate solution y™ of 
(A2) has the form, 

y n+1 =y n +u; n -5y n , (A3) 
where ui n is a parameter, Lu n < 1, and the correction 5y n is the solution of the matrix equation 

J™ • 5y n = -F n . (AA) 

In (A4) J is the Jacobian matrix, J; TO = dFi/dy m . The parameter u should be chosen to optimize the 
convergency of the iteration process to a solution. We use the following form of ui, 

" = 7 — aT> ( A5 ) 

max(i], A) 

where the parameter r\ = 0.03 and A is the average relative correction, 

1 K 
A=-Y 



Syk 



K 
k=i 



Vk 



In the case of equations (17), (18) we have K = 2(1 — 1) functional relations involved 21 variables, 
Vi and (c s )i- Two variables, vj and (c s )j, must be fixed as boundary conditions, when the number of 
independent variables equals to the number of equations, and system (A2) can be solved. In the case of 
equations (26)-(28) we have K = 3(1 — 1) and three boundary values, vj, (c s )j and Oj, which are fixed to 
provide a consistent solution of (A2). 

The presence of the singular points (r s ) in which coincides with n introduces additional complications 
to our method. To satisfy to two regularity conditions, 



D = and N = 0, 
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at n we modify the iteration procedure by the following way. We add the equation D = to the set of 
equations (Al) together with a new independent variable li n . Applying the Newton- Raphson iteration 
scheme we obtain a solution in which D = at r\, but N ^ in general. To satisfy the condition N — at 
the point i = 1 the appropriate choice of r\ must be done. We find the correct value of r\ using the bisection 
method in which n is changed by displacing all the grid. We apply the Newton-Raphson iterations (A4) 
for each new grid location. 

To correctly approximate the differential equations in the first grid interval (r\,r 2 ), where r\ — (r s ) in , 
we expand solution at r\ , 

v(r) =vi+v' s (r -n), (A7) 

c 2 s (r) = (c 2 s ) 1 + (cly s (r-n). (A8) 

The procedure of calculation of the coefficients v' s and (<? s )' s are given in Appendix B. Using (A7) and (A8) 
we fix the values at the second grid point as follows, 

v 2 =v 1 +v' s (r 2 -n), (A9) 

(c 2 s ) 2 = (cl) 1 + (cly s (r 2 -r 1 ). (AW) 

Relations (A9) and (A10) are used instead of the correspondent difference equations in (Al). The modified 
system of difference equations avoids numerical instabilities connected with the presence of the inner 
singular point and allows us to obtain a solution crossed continuously this point. 

We have found that described numerical procedure becomes unstable in the case of equations (26)-(28) 
at large values of a and m. The numerical instability arises because of influence of equation (26) and 
results in small-scale oscillations of all quantities. The numerical instability of similar type was found 
by Beloborodov (1998), who suppressed the oscillations on the level of ~ 10~ 3 of relative amplitude by 
applying a smoothing procedure to Ct. Unfortunately, such a smoothing procedure can not be included into 
our method because of additional coupling of equations (Al) with the regularity conditions (A6). 



B. Expansion at singular point 

We consider first the case of equations (26)-(28). We expand the numerator N 2 and denominator D 2 
at the singular point r s , as follows 

ft«-(^ + ^ + ^)(r-r.). m 

In (Bl) and (B2) the partial derivatives are taken at r = r s , v — v s , c 2 s — (Cg) s and fl = Cl s , and we use 
notations v s — v(r s ), (Cg) s = cj.(r s ) and fl s — Ct(r s ). The 'prime' mark means the radial derivative of the 
correspondent quantity in (Bl) and (B2). In (B2) we take into account that dD 2 /dfl — 0. We denote for 
convenience, 

< _ Mil - 

V s (cj) s il s 
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The value of Xs can be denned with help of (26). Substituting (Bl) and (B2) into equations (27) and (28) 
we finally obtain a quadratic equation with respect to £ s : 



dv 



del 



dD 2 
dr 



■2(c*) f 



dD 2 
del 



B 



dN 2 
del 



A 



dN 2 
dv 



is- 



dN 2 n dN 2 „. 2 , 



dN 2 
: del 



B 



(S3) 



where we denote 



A = l-Ml 



and 



2 r s 



+ 



n?r?- 



{r a -2f 



Having £ s from (B3) one obtain rj s = A£ s + B. The found values of £ s , ij s and Xs determine derivatives v' s , 
(Cg)g and fig which we use in the numerical procedure discussed in Appendix A. Equation (B3) has two 
roots. Which root should be used is defined by convergency of the iterations. We have found that only 
one of the two roots corresponds to the convergent solution. The partial derivatives of 7V 2 and D 2 used 
in (B3) can be calculated analytically or using numerical differentiation from (29) and (30). The latter 
method is simpler, but less accurate than the analytic one. The analytic derivation requires some efforts 
and produces quite long expressions which we do not present here. In the numerical procedure we have 
used the analytical expressions for derivatives, but in addition, we have checked them using estimates from 
the numerical differentiation. 



To judge the type of critical point we follow the procedure described by Kato et al. (1998). We 
introduce a new variable r defined by 



dr 



dr 
~r~D~ 2 



From equation (26)-(28) one obtain 



^= X VD 2 , ^=vN 2 , ^ = 2cl{AN 2 + BD 2 ). 
dr dr dr 



(B4) 



(B5) 



All of the variables are now expanded around those at r s , 



Ar, 



v s + Av, 



( C 2 s ) s + Ac 2 s , n = n s + An. 



(56) 



Substituting (B6) into (B4) and (B5) and retaining only linear terms, one obtain a system of linear 
differential equations with respect to Ar, Aw, Ac 2 and AO. Assuming that these quantities depend on r in 
the form exp(Ar), one obtain the following characteristic equation which determines the eigenvalues A: 



A — A 



dD 2 8N 2 . 2 . 
r„— +v a — + 2{c 2 s ) s 



2Ar s {c 2 s ) 



dr dv 
dN 2 8D 2 



A dN 2+B dD 2 



del 



8D 2 dN 2 + 8D 1 x ln dN 2 



dr dc 2 dr dc 2 dc 2 r, 
(dD 2 dN 2 dN 2 dD 2 



r,v. 



2Bv s (c 2 s)s v — — 
dD 2 dN 2 dN 2 dD 2 



dr dv 



dr dv 



dv dc 2 s 
dD 2 



dv 



dn 



dN 2 
dfl 



0. 



(57) 
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Again, all the partial derivatives in (B7) are taken at r = r s , v = v s , c? s = (c^) s and il = il s . If two solutions 
of quadratic equation (B7) are real and have different signs, A1A2 < 0, then the singular point is a saddle. 
The point is a node if two real roots of (B7) have identical signs, AiA 2 > 0. Complex conjugate roots of 
(B7) corresponds to the spiral singular point. 

In the case of equations (17), (18) equations (B3) and (B7) must be modified by substituting N\ and 
Di instead of N 2 and D 2 , and assuming dN\/dQ, = 0. The notation for B must be also changed to 
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Fig. 1. — Schematic illustration of behaviour of the integral curves in the vicinity of a global solution 
for transonic accretion disks. Two singular points exist at r = (r s )i n and (r s ) out in the case of viscosity 
prescription (1). The outer singular point at {r s ) out is always of a saddle type. The inner singular point 
at {r s ) in locates close to the last stable black hole orbit at 6r s and changes its type from a saddle to nodal 
one with increase of a as shown on panels (a) and (b) , respectively. Only the separatrix (thick line) which 
passes through both singular points represents the global solution. The nodal type inner singular point does 
not mean the existence of multiple physical solutions. 



(a) 



McVL E< 



(b) 




Fig. 2. — The specific angular momentum li n [panels (a) and (b)] and the position of the inner singular 
points [panels (c) and (d)] as a function of the mass accretion rate M for different viscosity parameters 
a = 0.01 (squares), 0.1 (circles) and 0.5 (triangles). Panels (a) and (c) correspond to viscosity prescriptions 
(1) and panels (b) and (d) correspond to viscosity prescriptions (2). The solid dots represent models with 
the saddle-type inner singular points, whereas the empty dots correspond to the nodal-type ones. 
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Fig. 3. — Location of the outer singular points as a function of the mass accretion rate M in the case 
of viscosity prescription (1). Models with a = 0.01, 0.1 and 0.5 arc shown. The points are always of a 
saddle-type. See caption to Fig. 2 for details of notations. 
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Fig. 4. — The specific angular momentum distribution with respect to radius in the innermost region of the 
disk with m = 160 and viscosity prescription (1). The long-dashed line corresponds to the Keplerian angular 
momentum. The short-dashed, solid and dotted lines correspond to a = 0.01, 0.1 and 0.5, respectively. The 
dots on the curves show the position of the inner singular points and the correspondent angular momentum. 
In the case of a = 0.01 and 0.1 the singular points are of a saddle-type, and in the case of a = 0.5 the point 
is of a nodal-type. 
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Fig. 5.— Ratio of the gas pressure to the total pressure at the inner singular point, (3 S , as a function of the 
mass accretion rate M. Panels (a) and (b) correspond to viscosity prescriptions (1) and (2), respectively. 
Models with a = 0.01, 0.1 and 0.5 are shown. See caption to Fig. 2 for details of notations. 
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Fig. 6. — Luminosities of accretion disks at different values of M and a in the case of viscosity prescription 
(1). The long-dashed line corresponds to the standard model with the radiative efficiency 1 / 16. Deviation 
of luminosities from one for the standard model indicates importance of advection. The short-dashed, solid 
and dotted lines correspond to a = 0.01, 0.1 and 0.5, respectively. 



